#! /usr/bin/env python
#
"""Classes for computing graph homology.
"""
__docformat__ = 'reStructuredText'

from native cimport dict, list

## logging subsystem

import itertools
import logging


## application-local imports

from aggregate cimport AggregateList
from combinatorics cimport bernoulli, factorial, minus_one_exp, Permutation, PermutationIterator, InplacePermutationIterator
from homology cimport *
from iterators cimport IndexedIterator
from rg cimport BoundaryCycle, Fatgraph, Vertex, MgnGraphsIterator
from rational cimport Rational
from runtime import runtime
from simplematrix import SimpleMatrix


NullMatrix = SimpleMatrix(0,0)


cpdef MgnChainComplex FatgraphComplex(int g, int n):
    """Return the fatgraph complex for given genus `g` and number of
    boundary components `n`.

    This is a factory method returning a `homology.ChainComplex`
    instance, populated with the correct vector spaces and
    differentials to compute the graph homology of the space
    `M_{g,n}`.
    """

    ## Minimum number of edges is attained when there's only one
    ## vertex; so, by Euler's formula `V - E + n = 2 - 2*g`, we get:
    ## `E = 2*g + n - 1`.
    cdef int min_edges = 2*g + n - 1
    logging.debug("  Minimum number of edges: %d", min_edges)
    
    ## Maximum number of edges is reached in graphs with all vertices
    ## tri-valent, so, combining Euler's formula with `3*V = 2*E`, we
    ## get: `E = 6*g + 3*n - 6`.  These are also graphs corresponding
    ## to top-dimensional cells.
    cdef int top_dimension = 6*g + 3*n - 6
    logging.debug("  Maximum number of edges: %d", top_dimension)

    # build chain complex
    cdef MgnChainComplex C
    C = MgnChainComplex(top_dimension)

    # gather graphs
    cdef chi = Rational(0)
    cdef int grade
    cdef Fatgraph graph
    for graph in MgnGraphsIterator(g,n):
        grade = graph.num_edges - 1
        pool = NumberedFatgraphPool(graph)
        # compute orbifold Euler characteristics (needs to include *all* graphs)
        chi += Rational(minus_one_exp(grade-min_edges)*len(pool), pool.num_automorphisms)
        if not pool.is_orientable:
            continue # discard non-orientable graphs
        C.module[grade].aggregate1(pool)
    C.orbifold_euler_characteristics = chi
    return C


cdef class MgnChainComplex(ChainComplex):
    """A specialized `ChainComplex`.
    """
    def __init__(self, int length):
        cdef int i
        ChainComplex.__init__(self, length)
        for i in xrange(length):
            self.module[i] = AggregateList()

    cpdef DifferentialComplex compute_boundary_operators(self):
        #: Matrix form of boundary operators; the `i`-th differential
        #: `D[i]` is `dim C[i-1]` rows (range) by `dim C[i]` columns
        #: (domain).
        cdef bint swap_rows_and_cols
        cdef int i
        cdef int j
        cdef int j0
        cdef int k
        cdef int k0
        cdef int l
        cdef int p
        cdef int pp
        cdef int q
        cdef int qq
        cdef int s

        cdef DifferentialComplex D
        cdef list m
        cdef NumberedFatgraphPool pool1
        cdef NumberedFatgraphPool pool2
        cdef object d

        cdef object checkpoint
        
        m = self.module # micro-optimization
        D = DifferentialComplex()
        D.append(NullMatrix, 0, len(m[0]))
        for i in xrange(1, len(self)):
            p = len(m[i-1]) # == dim C[i-1]
            q = len(m[i])   # == dim C[i]

            # for getting the most from the rank computing algorithms,
            # the second dimension should always be *larger* (because
            # it maps to the innermost loop)
            pp = p
            qq = q
            if p > q:
                # swap `p` and `q`
                p, q = q, p
                # keep track of it
                swap_rows_and_cols = True
            else:
                swap_rows_and_cols = False
                
            # maybe load `D[i]` from persistent storage
            checkpoint = runtime.options.checkpoint_dir + '/' + "M%d,%d-D%d.txt" % (runtime.g, runtime.n, i)
            if p>0 and q>0 and runtime.options.restart:
                d = SimpleMatrix(p, q)
                if d.load(checkpoint):
                    D.append(d, pp, qq)
                    logging.info("  Loaded %dx%d matrix D[%d] from file '%s'", pp, qq, i, checkpoint)
                    continue # with next `i`

            # compute `D[i]`
            d = SimpleMatrix(p, q)
            j0 = 0
            for pool1 in m[i].iterblocks():
                k0 = 0
                for pool2 in m[i-1].iterblocks():
                    for l in xrange(pool1.graph.num_edges):
                        if pool1.graph.is_loop(l):
                            continue # with next `l`
                        for (j, k, s) in NumberedFatgraphPool.facets(pool1, l, pool2):
                            assert k < len(pool2)
                            assert j < len(pool1)
                            assert k+k0 < pp
                            assert j+j0 < qq
                            if swap_rows_and_cols:
                                d.addToEntry(j+j0, k+k0, s)
                            else:
                                d.addToEntry(k+k0, j+j0, s)
                    k0 += len(pool2)
                    pool2.graph._cache_isomorphisms.clear()
                j0 += len(pool1)
                pool1.graph._cache_isomorphisms.clear()
            d.save(checkpoint)
            D.append(d, pp, qq)
            logging.info("  Computed %dx%d matrix D[%d]", pp, qq, i)
            
        return D



cdef class NumberedFatgraph(Fatgraph):
    """A `Fatgraph` decorated with a numbering of the boundary components.

    A numbered fatgraph is constructed from a `Fatgraph` instance
    (called the *underlying graph*) and a numbering (that is, a
    bijective map assigning an integer to each boundary components of
    the underlying graph).

    Examples::

      >>> ug = Fatgraph([Vertex([1, 0, 1, 0])])  # underlying graph
      >>> ng = NumberedFatgraph(ug, \
                 numbering=[(BoundaryCycle([(0,3,0), (0,2,3), (0,1,2), (0,0,1)]), 0)])

    The `numbering` argument must be a valid `dict` initializer.  The
    initializer can be:
      - either a sequence of tuples `(bcy, n)`, where each `n` is
        a non-negative integer, and each `bcy` is a
        `BoundaryCycle` instance,
      - or a `dict` instance mapping `BoundaryCycle` instances to
        `int`s.
    In either case, an assertion is raised if:
      - the number of pairs in the initializer does not match
        the number of boundary cycles;
      - the set of integer keys is not `[0 .. n]`;
      - there are duplicate boundary cycles or integers
        in the initializer;

    The `numbering` attribute is set to a dictionary mapping the
    boundary cycle `bcy` to the integer `n`::

      >>> ug0 = Fatgraph([Vertex([1,2,0]), Vertex([1,0,2])])
      >>> bc = ug0.boundary_cycles  # three b.c.'s
      >>> ng0 = NumberedFatgraph(ug0, [ (bcy,n) for (n,bcy) in enumerate(bc)])
      >>> ng0.numbering == {
      ...    BoundaryCycle([(0,0,1), (1,2,0)]): 0, 
      ...    BoundaryCycle([(0,1,2), (1,1,2)]): 1, 
      ...    BoundaryCycle([(0,2,0), (1,0,1)]): 2,
      ... }
      True

    Since `NumberedFatgraphs` are just decorated `Fatgraphs`, they
    only differ in the way two `NumberedFatgraph` instances are deemed
    isomorphic:
      - they must have isomorphic underlying graphs;
      - the numberings must match under the isomorphism map.
    For example::
    
      >>> ug = Fatgraph([Vertex([1, 1, 0, 0])])
      >>> ng1 = NumberedFatgraph(ug, numbering=[(BoundaryCycle([(0,3,0), (0,1,2)]), 0), \
                                                (BoundaryCycle([(0,0,1)]), 1),          \
                                                (BoundaryCycle([(0,2,3)]), 2)])
      >>> ng2 = NumberedFatgraph(ug, numbering=[(BoundaryCycle([(0,3,0), (0,1,2)]), 1), \
                                                (BoundaryCycle([(0,0,1)]), 0),          \
                                                (BoundaryCycle([(0,2,3)]), 2)])
      >>> ng1 == ng2
      False

      Fatgraph instances equipped with a numbering are compared as
      numbered graphs (that is, the isomorphism should transform the
      numbering on the source graph onto the numbering of the
      destination)::

        >>> NumberedFatgraph.__eq__(
        ...     NumberedFatgraph(Fatgraph([Vertex([2,0,1]), Vertex([2,1,0])]), 
        ...                      numbering=[(BoundaryCycle([(0,2,0), (1,0,1)]), 0), 
        ...                                 (BoundaryCycle([(0,0,1), (1,2,0)]), 1), 
        ...                                 (BoundaryCycle([(0,1,2), (1,1,2)]), 2) ] ),
        ...     NumberedFatgraph(Fatgraph([Vertex([2,0,1]), Vertex([2,1,0])]), 
        ...                      numbering=[(BoundaryCycle([(0,2,0), (1,0,1)]), 0), 
        ...                                 (BoundaryCycle([(0,0,1), (1,2,0)]), 2), 
        ...                                 (BoundaryCycle([(0,1,2), (1,1,2)]), 1) ] ) )
        True

        >>> NumberedFatgraph.__eq__(
        ...     NumberedFatgraph(Fatgraph([Vertex([1, 0, 0, 2, 2, 1])]), 
        ...                      numbering=[(BoundaryCycle([(0,5,0)]), 0), 
        ...                                 (BoundaryCycle([(0,0,1), (0,2,3), (0,4,5)]), 1), 
        ...                                 (BoundaryCycle([(0,1,2)]), 3), 
        ...                                 (BoundaryCycle([(0,3,4)]), 2) ]),
        ...     NumberedFatgraph( 
        ...                      Fatgraph([Vertex([2, 2, 1, 1, 0, 0])]), 
        ...                      numbering=[(BoundaryCycle([(0,2,3)]), 0), 
        ...                                 (BoundaryCycle([(0,3,4), (0,5,0), (0,1,2)]), 3), 
        ...                                 (BoundaryCycle([(0,4,5)]), 1), 
        ...                                 (BoundaryCycle([(0,0,1)]), 2) ]) )
        False

        >>> NumberedFatgraph.__eq__(
        ...     NumberedFatgraph(Fatgraph([Vertex([1, 0, 0, 2, 2, 1])]), 
        ...                      numbering=[(BoundaryCycle([(0,5,0)]), 0), 
        ...                                 (BoundaryCycle([(0,0,1), (0,2,3), (0,4,5)]), 1),
        ...                                 (BoundaryCycle([(0,1,2)]), 3), 
        ...                                 (BoundaryCycle([(0,3,4)]), 2) ]),
        ...     NumberedFatgraph(Fatgraph([Vertex([2, 2, 1, 1, 0, 0])]), 
        ...                      numbering=[(BoundaryCycle([(0,2,3)]), 3), 
        ...                                 (BoundaryCycle([(0,3,4), (0,5,0), (0,1,2)]), 2), 
        ...                                 (BoundaryCycle([(0,4,5)]), 0), 
        ...                                 (BoundaryCycle([(0,0,1)]), 1) ]) )
        False

        >>> NumberedFatgraph.__eq__(
        ...     NumberedFatgraph(Fatgraph([Vertex([3, 2, 2, 0, 1]), Vertex([3, 1, 0])]), 
        ...                      numbering=[(BoundaryCycle([(0,4,0), (1,0,1)]), 0),
        ...                                 (BoundaryCycle([(0,0,1), (0,2,3), (1,2,0)]), 1),
        ...                                 (BoundaryCycle([(0,1,2)]), 2),
        ...                                 (BoundaryCycle([(0,3,4), (1,1,2)]), 3) ]),
        ...     NumberedFatgraph(Fatgraph([Vertex([2, 3, 1]), Vertex([2, 1, 3, 0, 0])]), 
        ...                      numbering=[(BoundaryCycle([(0,2,0), (1,0,1)]), 3),
        ...                                 (BoundaryCycle([(0,0,1), (1,2,3), (1,4,0)]), 1),
        ...                                 (BoundaryCycle([(0,1,2), (1,1,2)]), 0) ,
        ...                                 (BoundaryCycle([(1,3,4)]), 2) ]) )
        True

        >>> NumberedFatgraph.__eq__(
        ...     NumberedFatgraph(Fatgraph([Vertex([0, 1, 2, 0, 2, 1])]),
        ...                      numbering=[(BoundaryCycle([(0,1,2), (0,4,5)]), 0),
        ...                                 (BoundaryCycle([(0,5,0), (0,3,4), (0,2,3), (0,0,1)]), 1)]),
        ...     NumberedFatgraph(Fatgraph([Vertex([0, 1, 2, 0, 2, 1])]),
        ...                      numbering=[(BoundaryCycle([(0,1,2), (0,4,5)]), 1),
        ...                                 (BoundaryCycle([(0,5,0), (0,3,4), (0,2,3), (0,0,1)]), 0)]) )
        False
      """

    def __init__(self, underlying, numbering):
        Fatgraph.__init__(self, underlying)
        self.underlying = underlying
        assert len(numbering) == self.num_boundary_cycles
        self.numbering = dict(numbering)


    def __repr__(self):
        """Output a printed representation, such that `eval(repr(x)) == x`.

        The `numbering` attribute of the `NumberedFatgraph` differs from
        the standard Python printing of dictionaries, in that its printed
        form is sorted by values (i.e., by boundary component index)
        to make doctests more stable.
        """
        cdef list parts
        parts = list()
        for (k, v) in sorted(self.numbering.iteritems()):
            parts.append("%s: %s" % (repr(k), repr(v)))
        canonical = "{" + str.join(", ", parts) + "}"
        return "NumberedFatgraph(%s, numbering=%s)" \
               % (repr(self.underlying), canonical)
    

    #@ocache_weakref
    cdef NumberedFatgraph contract_nfg(self, int edgeno):
        """Return a new `NumberedFatgraph` instance, obtained by
        contracting the specified edge.

        Examples::

          >>> BoundaryCycle = BoundaryCycle
          >>> g0 = NumberedFatgraph(Fatgraph([Vertex([1, 2, 1]), Vertex([2, 0, 0])]),
          ...                       numbering={BoundaryCycle([(0, 1, 2), (1, 2, 0),
          ...                                                 (0, 0, 1), (1, 0, 1)]): 0,
          ...                                  BoundaryCycle([(1, 1, 2)]): 1,
          ...                                  BoundaryCycle([(0, 2, 0)]): 2})
          >>> g0.contract(2)
          NumberedFatgraph(Fatgraph([Vertex([1, 1, 0, 0])]),
                           numbering={BoundaryCycle([(0, 3, 0), (0, 1, 2)]): 0,
                                      BoundaryCycle([(0, 2, 3)]): 1,
                                      BoundaryCycle([(0, 0, 1)]): 2})

          >>> g1 = NumberedFatgraph(Fatgraph([Vertex([1, 0, 2]), Vertex([2, 1, 5]),
          ...                                 Vertex([0, 4, 3]), Vertex([4, 5, 3])]),
          ...                       numbering={BoundaryCycle([(0, 2, 0), (0, 1, 2), (0, 0, 1),
          ...                                                 (1, 1, 2), (2, 0, 1), (3, 0, 1),
          ...                                                 (1, 2, 0), (2, 2, 0), (1, 0, 1),
          ...                                                                       (3, 1, 2)]): 0,
          ...                                  BoundaryCycle([(3, 2, 0), (2, 1, 2)]): 1})
          >>> g1.contract(0)
          NumberedFatgraph(Fatgraph([Vertex([1, 0, 3, 2]), Vertex([1, 0, 4]), Vertex([3, 4, 2])]),
                           numbering={BoundaryCycle([(2, 1, 2), (0, 1, 2), (0, 0, 1), (1, 1, 2),
                                                     (0, 3, 0), (2, 0, 1), (1, 2, 0), (1, 0, 1)]): 0,
                                      BoundaryCycle([(0, 2, 3), (2, 2, 0)]): 1})

        """
        # check that we are not contracting a loop or an external edge
        assert not self.is_loop(edgeno), \
               "NumberedFatgraph.contract: cannot contract a loop."
        assert (edgeno >= 0) and (edgeno < self.num_edges), \
               "NumberedFatgraph.contract: invalid edge number (%d):"\
               " must be in range 0..%d" \
               % (edgeno, self.num_edges)

        # store endpoints of the edge-to-be-contracted
        cdef int v1
        cdef int v2
        cdef int pos1
        cdef int pos2
        cdef int l1
        cdef int l2
        cdef int n
        cdef int i1
        cdef int i2
        cdef dict new_numbering
        cdef BoundaryCycle bcy
        cdef list new_cy
        cdef Fatgraph contracted

        (v1, v2) = self.endpoints_v[edgeno]
        (pos1, pos2) = self.endpoints_i[edgeno]
        if v1 > v2:
            # swap endpoints so that `v1 < v2`
            v1, v2 = v2, v1
            pos1, pos2 = pos2, pos1
        l1 = (len(self.vertices[v1]))
        l2 = (len(self.vertices[v2]))
        
        # transform corners according to contraction; see
        # `Fatgraph.contract()` for an explanation of how the
        # underlying graph is altered during contraction.
        contracted = self.underlying.contract_fg(edgeno)
        new_numbering = dict()
        for (bcy, n) in self.numbering.iteritems():
            new_cy = bcy.contract_bcy((v1,pos1), (v2,pos2), contracted)
            new_numbering[Fatgraph.BoundaryCycle(new_cy)] = n
        return NumberedFatgraph(contracted, numbering=new_numbering)

        
    cpdef list isomorphisms(self, other):
        """Return list of isomorphisms from `self` to `other`.

        See `Fatgraph.isomrphisms` for a discussion of the
        representation of isomorphisms and example usage.
        """
        # try cached result first
        key = other.cache_id()
        try:
            return self._cache_isomorphisms[key]
        except KeyError:
            pass

        cdef list result
        cdef dict pe
        cdef dict pv
        cdef list rot
        cdef BoundaryCycle bc1
        cdef BoundaryCycle bc2
        cdef bint pe_does_not_preserve_bc
        
        result = list()
        for (pv, rot, pe) in Fatgraph.isomorphisms(self.underlying, other.underlying):
            pe_does_not_preserve_bc = False
            for bc1 in self.underlying.boundary_cycles:
                bc2 = bc1.transform((pv, rot, pe), other.underlying)
                # there are cases (see examples in the
                # `Fatgraph.__eq__` docstring, in which the
                # above algorithm may find a valid
                # mapping, changing from `g1` to an
                # *alternate* representation of `g2` -
                # these should fail as they don't preserve
                # the boundary cycles, so we catch them
                # here.
                if (bc2 not in other.numbering) \
                       or (self.numbering[bc1] != other.numbering[bc2]):
                    pe_does_not_preserve_bc = True
                    break
            if pe_does_not_preserve_bc:
                continue # to next underlying graph isomorphism
            result.append((pv, rot, pe))

        # at last
        self._cache_isomorphisms[key] = result
        return result



# ---- auxiliary functions for NumberedFatgraphPool

cdef inline bint unseen(list candidate, list P, list already):
    """Return `False` iff any of the images of `candidate` by an
    element of group `P` is contained in set `already`.
    """
    for p in P:
        if p.rearranged(candidate) in already:
            return False
    return True

cdef inline Permutation push_fwd(f, Fatgraph g1, Fatgraph g2):
    """Return a `Permutation` instance corresponding to the
    map induced by `f` on the boundary cycles, or `None`
    if no such map exists.

    Indeed, there are cases (see examples in the
    `Fatgraph.__eq__` docstring), in which the
    `Fatgraph.isomorphisms` algorithm may find a valid
    mapping, changing from `g1` to an *alternate*
    representation of `g2` - these should fail as they don't
    preserve the boundary cycles.
    """
    cdef Permutation result
    cdef int src
    cdef BoundaryCycle bc1
    cdef BoundaryCycle bc2
    
    result = Permutation() #: map induced by `f` on boundary cycles
    for (src, bc1) in enumerate(g1.boundary_cycles):
        bc2 = bc1.transform(f, g2)
        try:
            result[src] = g2.boundary_cycles.index(bc2)
        except KeyError:
            return None
    return result

cdef inline list pull_back(Permutation f0, Permutation f1, list nb):
    return f1.rearranged(f0.inverse().rearranged(nb))

cdef inline compute_nb_map(Permutation f0, Permutation f1,
                    NumberedFatgraphPool src, NumberedFatgraphPool dst):
    """
    Return a pair of `dict` instances:
      - the first maps (indices of ) numberings in
        `NumberedFatgraphPool` instance `src` to corresponding
        (indices of) numberings in `NumberedFatgraphPool`
        instance `dst`.

      - the second maps the same indices into the
        automorphism `a` of `dst.graph` such that::

          src.numberings[i] = pull_back(<permutation induced by a applied to> dst.numberings[i])

    Raises `KeyError` if no such mapping could be found.

    Indices must be `int`s here, as `list`s cannot be used as
    dictionary keys.
    """
    cdef dict nb_map
    cdef dict a_map
    cdef int i
    cdef int i_
    cdef list nb
    cdef object k
    
    nb_map = dict()
    a_map = dict()
    for (i, nb) in enumerate(src.numberings):
        (i_, k) = dst._index(pull_back(f0, f1, nb))
        nb_map[i] = i_
        a_map[i] = k
    return (nb_map, a_map)

# ---- end of auxiliary functions for NumberedFatgraphPool

cdef class NumberedFatgraphPool(object):
    """An immutable virtual collection of `NumberedFatgraph`s.
    Items are all distinct (up to isomorphism) decorations of a
    `Fatgraph` instance `graph` with a numbering of the boundary
    cycles.

    Implements object lookup by index and iteration over the whole list.

    Examples::
    
      >>> ug1 = Fatgraph([Vertex([2,0,0]), Vertex([2,1,1])])
      >>> p = NumberedFatgraphPool(ug1)
      >>> len(p)
      3
      >>> for g in p: print g
      NumberedFatgraph(Fatgraph([Vertex([2, 0, 0]), Vertex([2, 1, 1])]),
                       numbering={BoundaryCycle([(0, 2, 0), (1, 2, 0), (0, 0, 1), (1, 0, 1)]): 0,
                                  BoundaryCycle([(1, 1, 2)]): 1,
                                  BoundaryCycle([(0, 1, 2)]): 2})
      NumberedFatgraph(Fatgraph([Vertex([2, 0, 0]), Vertex([2, 1, 1])]),
                       numbering={BoundaryCycle([(0, 1, 2)]): 0,
                                  BoundaryCycle([(0, 2, 0), (1, 2, 0), (0, 0, 1), (1, 0, 1)]): 1,
                                  BoundaryCycle([(1, 1, 2)]): 2})
      NumberedFatgraph(Fatgraph([Vertex([2, 0, 0]), Vertex([2, 1, 1])]),
                       numbering={BoundaryCycle([(0, 1, 2)]): 0,
                                  BoundaryCycle([(1, 1, 2)]): 1,
                                  BoundaryCycle([(0, 2, 0), (1, 2, 0), (0, 0, 1), (1, 0, 1)]): 2})
       
    Note that, when only one numbering out of many possible ones is
    returned because of isomorphism, the returned numbering may not be
    the trivial one (it is infact the first permutation of 0..n
    returned by `InplacePermutationIterator`)::
      
      >>> ug2 = Fatgraph([Vertex([2,1,0]), Vertex([2,0,1])])
      >>> for g in NumberedFatgraphPool(ug2): print g
      NumberedFatgraph(Fatgraph([Vertex([2, 1, 0]), Vertex([2, 0, 1])]),
                        numbering={BoundaryCycle([(1, 2, 0), (0, 0, 1)]): 0,
                                   BoundaryCycle([(0, 2, 0), (1, 0, 1)]): 1,
                                   BoundaryCycle([(0, 1, 2), (1, 1, 2)]): 2})

    When the graph has only one boundary component, there is only one
    possible numbering, which is actually returned::
    
      >>> ug3 = Fatgraph([Vertex([1,0,1,0])])
      >>> for g in NumberedFatgraphPool(ug3): print g
      NumberedFatgraph(Fatgraph([Vertex([1, 0, 1, 0])]),
                        numbering={BoundaryCycle([(0, 3, 0), (0, 2, 3),
                                                  (0, 1, 2), (0, 0, 1)]): 0})

    Index lookup returns a `NumberedFatgraph` instance, produced
    on-the-fly::

      >>> ug4 = Fatgraph([Vertex([0, 1, 0, 1, 2, 2])])
      >>> pool = NumberedFatgraphPool(ug4)
      >>> pool[0]
      NumberedFatgraph(Fatgraph([Vertex([0, 1, 0, 1, 2, 2])]),
                       numbering={BoundaryCycle([(0, 4, 5)]): 0,
                                  BoundaryCycle([(0, 2, 3), (0, 3, 4),
                                                 (0, 1, 2), (0, 0, 1), (0, 5, 0)]): 1})
    """
    def __init__(self, graph):
        cdef list bc
        cdef int n
        cdef bint orientable
        cdef list P
        cdef list P_
        cdef list automorphisms
        cdef int src
        cdef BoundaryCycle dst_cy
        cdef int dst
        cdef Permutation p        
        cdef list numberings
        cdef list candidate
        
        bc = graph.boundary_cycles
        n = len(bc) # == graph.num_boundary_cycles
        orientable = True

        if n>1:
            ## Find out which automorphisms permute the boundary cycles among
            ## themselves.
            P = list()  #: permutation of boundary cycles induced by `a \in Aut(G)`
            P_ = list() #: corresponding graph automorphisms: `P[i]` is induced by `P_[i]`
            automorphisms = list() #: `NumberedFatgraph` automorphisms
            for a in graph.automorphisms():
                p = Permutation()
                for src in xrange(n):
                    dst_cy = bc[src].transform(a, graph)
                    try:
                        dst = bc.index(dst_cy)
                    except ValueError:
                        # `dst_cy` not in `bc`
                        break # continue with next `a`
                    p[src] = dst
                if len(p) != n: # not all `src` were mapped to a `dst`
                    continue # with next `a`
                if p.is_identity():
                    # `a` preserves the boundary cycles pointwise,
                    # so it induces an automorphism of the numbered graph
                    automorphisms.append(a)
                    if Fatgraph._cmp_orient(graph, graph, a) == -1:
                        orientable = False
                if (p not in P):
                    # `a` induces permutation `p` on the set `bc`
                    P.append(p)
                    P_.append(a)
            assert len(P) > 0 # XXX: should verify that `P` is a group!
        else: # n==1
            # there's no permutation of boundary cycles except the trivial one,
            # which is clearly induced by the trivial automorphism
            P = list([ Permutation([(0,0)]) ])
            P_ = list()
            # some Fatgraph automorphisms might map to a different representation
            # of the same graph; those need to be discarded
            automorphisms = list()
            for a in graph.automorphisms():
                dst_cy = bc[0].transform(a, graph)
                if dst_cy != bc[0]:
                    continue # with next `a`
                if Fatgraph._cmp_orient(graph, graph, a) == -1:
                    orientable = False
                if len(P_) == 0:
                    P_.append(a)
                automorphisms.append(a)

        ## There will be as many distinct numberings as there are cosets
        ## of `P` in `Sym(n)`.
        if len(P) > 1:
            numberings = list()
            for candidate in InplacePermutationIterator(range(n)):
                if unseen(candidate, P, numberings):
                    numberings.append(candidate[:])
        else:
            # if `P` is the one-element group, then all orbits are trivial
            #numberings = PermutationList(n)
            #numberings = list(PermutationIterator(range(n)))
            numberings = [ candidate[:] for candidate in InplacePermutationIterator(range(n)) ]

        # things to remember
        self.graph = graph
        self.is_orientable = orientable
        self.P = P
        self.P_ = P_
        self.numberings = numberings
        self.num_automorphisms = len(automorphisms)


    def __getitem__(self, pos):
        return NumberedFatgraph(self.graph, zip(self.graph.boundary_cycles, self.numberings[pos]))


    def __iter__(self):
        return IndexedIterator(self)
        

    def __len__(self):
        return len(self.numberings)


    def __repr__(self):
        if hasattr(self, 'graph'):
            return "NumberedFatgraphPool(%s)" % self.graph
        else:
            return object.__repr__(self)
    def __str__(self):
        return repr(self)
        

    cdef list facets(self, int edge, NumberedFatgraphPool other):
        """Iterate over facets obtained by contracting `edge` and
        projecting onto `other`.

        Each returned item is a triple `(j, k, s)`, where:
          - `j` is the index of a `NumberedFatgraph` in `self`;
          - `k` is the index of a `NumberedFatgraph` in `other`;
          - `s` is the sign by which `self[j].contract(edge)` projects onto `other[k]`.
        Only triples for which `s != 0` are returned.

        Examples::
        
          >>> p0 = NumberedFatgraphPool(Fatgraph([Vertex([1, 2, 0, 1, 0]), Vertex([3, 3, 2])]))
          >>> p1 = NumberedFatgraphPool(Fatgraph([Vertex([0, 1, 0, 1, 2, 2])]))
          >>> list(NumberedFatgraphPool.facets(p0, 2, p1))
          [(0, 0, 1), (1, 1, 1)]
        """
        cdef Fatgraph g0
        cdef Fatgraph g1
        cdef Fatgraph g2
        cdef Permutation f0_push_fwd
        cdef int i
        cdef int j
        cdef int k
        cdef int s
        cdef list nb
        cdef bint f1_found
        cdef Permutation f1_push_fwd
        cdef object a  # tuple
        cdef object e1 # tuple
        cdef object e2 # tuple
        cdef list result
        
        assert not self.graph.is_loop(edge)
        assert self.is_orientable
        assert other.is_orientable
        
        g0 = self.graph
        g1 = g0.contract_fg(edge)
        g2 = other.graph
        assert len(g1.boundary_cycles) == len(g2.boundary_cycles)
        
        ## 1. compute map induced on `g0.boundary_cycles` from the
        ##    graph map `f0` which contracts `edge`.
        (e1, e2) = g0.endpoints(edge)
        assert set(g1.boundary_cycles) == set([ bcy.contract_bcy(e1, e2, g1) for bcy in g0.boundary_cycles ]), \
               "NumberedFatgraphPool.facets():" \
               " Boundary cycles of contracted graph are not the same" \
               " as contracted boundary cycles of parent graph:" \
               " `%s` vs `%s`" % (g1.boundary_cycles, [ bcy.contract_bcy(e1, e2, g1)
                                                        for bcy in g0.boundary_cycles ])
        f0_push_fwd = Permutation(enumerate([
            g1.boundary_cycles.index(bcy.contract_bcy(e1, e2, g1))
            for bcy in g0.boundary_cycles
            ]))
        assert len(f0_push_fwd) == len(g1.boundary_cycles)
        assert len(f0_push_fwd) == g0.num_boundary_cycles

        ## 2. compute isomorphism map `f1` from `other.graph` to `self.graph.contract(edge)`;
        ##    if there is no such isomorphisms, then stop iteration.
        ##    The reason for having this map go in the opposite direction as `f0`
        ##    is to avoid inverting its push-forward.  (Thus saving a little complexity.)

        # if there is no isomorphism from `g2` to `g1`, the "for" loop
        # down here will not run, so we need a way to detect that a
        # valid isomorphism has been found.
        f1_found = False
        for f1 in Fatgraph.isomorphisms(g2, g1):
            f1_push_fwd = push_fwd(f1, g2, g1)
            if f1_push_fwd is None:
                continue # with next isomorphism
            # now check that every numbering of `g0` can be mapped to
            # a numbering on `g2` through the composite map `f1^(-1) * f0`
            try:
                nb_map, a_map = compute_nb_map(f0_push_fwd, f1_push_fwd, self, other)
            except KeyError:
                continue # with next isomorphism
            # if we got to this point, the `f1` is a graph
            # automorphism that induces a map on the boundary cycles
            f1_found = True
            break # correct `f1` found
        if not f1_found:
            # no isomorphism between `g1` and `g2`, stop here.
            return list()

        ## For every numbering `nb` on `g0`, compute the (index of)
        ## corresponding numbering on `g2` (under the composition map
        ## `f1^(-1) * f0`) and return the triple `(index of nb, index of
        ## push-forward, sign)`.
        result = list()
        for (j, nb) in enumerate(self.numberings):
            k = nb_map[j]
            a = a_map[j]
            ## there are three components to the sign `s`:
            ##   - the sign given by the ismorphism `f1`
            ##   - the sign of the automorphism of `g2` that transforms the
            ##     push-forward numbering into the chosen representative in the same orbit
            ##   - the alternating sign from the homology differential
            s = Fatgraph._cmp_orient(g2, g1, f1) \
                * Fatgraph._cmp_orient(g2, g2, a) \
                * minus_one_exp(g0.edge_numbering[edge])
            result.append((j, k, s))

        return result
    

    cdef _index(self, list numbering):
        """
        """
        cdef int i
        cdef Permutation p
        cdef int j
        
        for (i, p) in enumerate(self.P):
            try:
                j = self.numberings.index(p.rearranged(numbering))
                # once a `p` has matched, there's no reason to try others
                return (j, self.P_[i])
            except ValueError:
                pass
        assert False, \
               "%s._index(%s): No match found." % (self, numbering)
